/**
 * @defgroup   GLNI
 *
 * @brief	  This head file contains implementations of the Gauss Legendre Numerical Integral on 1D, 2D and 3D ranges.
 * 
 * This head file defines five classes, including the line, quadrilateral, triangle, cube and tetrahedron for the GLNI on 
 * corresponding 2D and 3D ranges. All classes provide integrals of customized functions over defined 
 * ranges though a callback interface up to 20 orders of the GLNI. 
 * All classes preform investigations of the parameters' validity by default. The user could define the GLNI_NO_CHECK 
 * macro to disable this functionality, which might speed up the integrals. Detailed documentations of different 
 * classes are given below. If you want to know more, please contact the author at zhangyiss@icloud.com.
 *
 * @author	 Dr. Yi Zhang
 * @date	   2021-03-14
 * @email	  zhangyiss@icloud.com
 */
#ifndef _GLNI_H
#define _GLNI_H

#include "exception"
#include "stdexcept"
#include "complex"
#include "vector"

#ifndef GLNI_ZERO
#define GLNI_ZERO 1e-10
#endif // GLNI_ZERO

namespace glni
{


double interpolate1d_linear(double x1, double x2, double y1, double y2, double x);

/**
 * @brief	Gauss Legendre Numerical Integration in the reference 1D segment
 * 
 * integral range for the reference element: [-1, 1]
 * -1         1
 * |----------|----> u
 * 
 * a class object must be initiated before use. This could be done by passing a desired node 
 * number of the GLNI to the constructor or by calling the initiate() function. To use 
 * a different order the GLNI, one must firstly calling the clear() function before the class 
 * object could be initiated again. To use the class object for integration, the recommended way 
 * is to use the predefined callback interface for customized functions. Once the functions are 
 * defined. It could be passed to the integral() function for evaluation. Also, one could use the 
 * node() and weight() functions to get the Gauss nodes' positions and weights, respectively. 
 * Note that the nodes' position are located wrt. the reference element. Proper scaling must be 
 * applied to map the nodes to a customized range.
 */
class line
{
public:
	line();
	line(unsigned int size);
	virtual ~line();

	void clear();
	int size(){return size_;}
	double node(unsigned int idx);
	double weight(unsigned int idx);

	typedef double (*func_1d_ptr)(double x, double x1, double x2, void *att, void *instance);

	typedef double (*func_1d_ptr2)(double x, const std::vector<double> &ix, const std::vector<double> &iv, void *att, void *instance);

	/**
	 * @brief      Integration over a 1D range
	 *
	 * @param[in]  low       The lower bound
	 * @param[in]  hig       The higher bound
	 * @param[in]  fx        The function to be integrated
	 * @param      att       Additional attributes need to be passed to the function
	 * @param      instance  The user object passed to the function
	 *
	 * @return     Integral value
	 */
	double integral(double low, double hig, func_1d_ptr fx, void *att = nullptr, void *instance = nullptr);

	/**
	 * @brief      Integration over a 1D range
	 * 
	 * @param[in] ix             x coordinates
	 * @param[in] iv             corresponding measurements
	 * @param[in] fx             The function to be integrated
	 * @param     att            Additional attributes need to be passed to the function
	 * @param     instance       The user object passed to the function
	 * @return    Integral value
	 */
	double integral(const std::vector<double> &ix, const std::vector<double> &iv, func_1d_ptr2 fx, void *att = nullptr, void *instance = nullptr);

	void initiate(unsigned int size);

protected:
	bool initialized_;
	int size_;
	double *nodes_;
	double *weights_;
};

/**
 * @brief	Gauss Legendre Numerical Integration in the 2D quadrilateral
 * 
 * integral range for the reference element: [-1, 1]x[-1, 1]
 * y
 * |
 * [-1,1]	 [1,1]
 * |----------|
 * |		  |
 * |		  |
 * |		  |
 * |----------|----> x
 * [-1,-1]   [1,-1]
 * 
 * a class object must be initiated before use. This could be done by passing a desired node 
 * number of the GLNI to the constructor or by calling the initiate() function. To use 
 * a different order the GLNI, one must firstly calling the clear() function before the class 
 * object could be initiated again. To use the class object for integration, the recommended way 
 * is to use the predefined callback interface for customized functions. Once the functions are 
 * defined. It could be passed to the integral() function for evaluation. Since specified transformations 
 * are required to calculate the integral in a 2D quadrilateral area. No direct access to the 
 * nodes and weights of the GLNI are provided here.
 */
class quadrilateral
{
public:
	quadrilateral();
	quadrilateral(unsigned int size);
	virtual ~quadrilateral();

	void clear();
	int size(){return g1d_.size();}

	typedef double (*func_quad_ptr)(double x, double y, double x1, double x2, double y1, double y2, void *att, void *instance);
	double integral(double xlow, double xhig, double ylow, double yhig, func_quad_ptr fx, void *att = nullptr, void *instance = nullptr);

	void initiate(unsigned int size);

protected:
	line g1d_;
	bool initialized_;
	int size_;
};

/**
 * @brief	Gauss Legendre Numerical Integration in the 2D triangle
 * 
 * integral range for the reference element: [0, 0]x[0, 1]x[0, 1]
 * v
 * |
 * [0, 1]
 * |\
 * |  \
 * |	\
 * |	  \
 * |		\
 * |		  \
 * |------------\---> u
 * [0, 0]     [1, 0]
 * 
 * a class object must be initiated before use. This could be done by passing a desired node 
 * number of the GLNI to the constructor or by calling the initiate() function. To use 
 * a different order the GLNI, one must firstly calling the clear() function before the class 
 * object could be initiated again. To use the class object for integration, the recommended way 
 * is to use the predefined callback interface for customized functions. Once the functions are 
 * defined. It could be passed to the integral() function for evaluation. Since specified transformations 
 * are required to calculate the integral in a 2D triangular area. No direct access to the 
 * nodes and weights of the GLNI are provided here.
 */
class triangle
{
public:
	triangle();
	triangle(unsigned int size);
	virtual ~triangle();

	void clear();
	int size(){return g1d_.size();}

	typedef double (*func_tri_ptr)(double x, double y, double x1, double x2, double x3, double y1, double y2, double y3, void *att, void *instance);
	double integral(double x1, double x2, double x3, double y1, double y2, double y3, func_tri_ptr fx, void *att = nullptr, void *instance = nullptr);

	typedef std::complex<double> (*func_tri_ptr2)(double x, double y, double x1, double x2, double x3, double y1, double y2, double y3, void *att, void *instance);
	std::complex<double> integral(double x1, double x2, double x3, double y1, double y2, double y3, func_tri_ptr2 fx, void *att = nullptr, void *instance = nullptr);

	void initiate(unsigned int size);

protected:
	line g1d_;
	bool initialized_;
	int size_;
};

/**
 * @brief	Gauss Legendre Numerical Integration in the reference 3D cube
 * 
 * integral range for the reference element: [-1, -1, -1] ~ [1, 1, 1]
 * z
 * |
 * |    ----------- [1,1,1]
 * |  /|         /|
 * | / |        / |
 * |/  |       /  |
 * -----------|   |
 * |   |      |   |
 * |  /-------|--/
 * | /        | /
 * |----------|--------> x
 * [-1,-1,-1]
 * 
 * a class object must be initiated before use. This could be done by passing a desired node 
 * number of the GLNI to the constructor or by calling the initiate() function. To use 
 * a different order the GLNI, one must firstly calling the clear() function before the class 
 * object could be initiated again. To use the class object for integration, the recommended way 
 * is to use the predefined callback interface for customized functions. Once the functions are 
 * defined. It could be passed to the integral() function for evaluation. Since specified transformations 
 * are required to calculate the integral in a 3D cubic space. No direct access to the 
 * nodes and weights of the GLNI are provided here.
 */
class cube
{
public:
	cube();
	cube(unsigned int size);
	virtual ~cube();

	void clear();
	int size(){return g1d_.size();}

	typedef double (*func_cube_ptr)(double x, double y, double z, 
		double x1, double x2, double y1, double y2, double z1, double z2, void *att, void *instance);
	double integral(double x1, double x2, double y1, double y2, 
		double z1, double z2, func_cube_ptr fx, void *att = nullptr, void *instance = nullptr);

	void initiate(unsigned int size);

protected:
	line g1d_;
	bool initialized_;
	int size_;
};

/**
 * @brief	Gauss Legendre Numerical Integration in the reference tetrahedron
 * 
 * integral range for the reference element: [0, 0, 0]x[1, 0, 0]x[0, 1, 0]x[0, 0, 1]
 * q
 * |
 * |  [0, 0, 1]
 * |\
 * |  \
 * |	\
 * |	  \
 * |		\
 * |		  \ [0, 1, 0]
 * |------------\---> v
 * | [0, 0, 0] /
 * |         /
 * |       /
 * |     /
 * |   /
 * | /
 * |/ [1, 0, 0]
 * |
 * |
 * u
 * 
 * a class object must be initiated before use. This could be done by passing a desired node 
 * number of the GLNI to the constructor or by calling the initiate() function. To use 
 * a different order the GLNI, one must firstly calling the clear() function before the class 
 * object could be initiated again. To use the class object for integration, the recommended way 
 * is to use the predefined callback interface for customized functions. Once the functions are 
 * defined. It could be passed to the integral() function for evaluation. Since specified transformations 
 * are required to calculate the integral in a 3D tetrahedral space. No direct access to the 
 * nodes and weights of the GLNI are provided here.
 */
class tetrahedron
{
public:
	tetrahedron();
	tetrahedron(unsigned int size);
	virtual ~tetrahedron();

	void clear();
	int size(){return g1d_.size();}

	typedef double (*func_tet_ptr)(double x, double y, double z, 
		double x1, double x2, double x3, double x4, 
		double y1, double y2, double y3, double y4,
		double z1, double z2, double z3, double z4, void *att, void *instance);
	double integral(double x1, double x2, double x3, double x4, 
		double y1, double y2, double y3, double y4, 
		double z1, double z2, double z3, double z4, 
		func_tet_ptr fx, void *att = nullptr, void *instance = nullptr);

	typedef std::complex<double> (*func_tet_ptr2)(double x, double y, double z, 
		double x1, double x2, double x3, double x4, 
		double y1, double y2, double y3, double y4,
		double z1, double z2, double z3, double z4, void *att, void *instance);
	std::complex<double> integral(double x1, double x2, double x3, double x4, 
		double y1, double y2, double y3, double y4, 
		double z1, double z2, double z3, double z4, 
		func_tet_ptr2 fx, void *att = nullptr, void *instance = nullptr);

	void initiate(unsigned int size);

protected:
	line g1d_;
	bool initialized_;
	int size_;
};

}

#endif //_GLNI_H